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Abstract 

Hadronic annihilation rate of l +_ heavy quarkonium is given to next-to-leading order in a s and 
leading order in v using a recently developed factorization formalism which is based on NRQCD. The 
result includes both the annihilation of P-wave color-singlet QQ component, and the annihilation of 
S-wave color-octet QQ component of the quarkonium. The notorious infrared divergences due to soft 
gluons, i.e., the Logarithms associated with the binding energy, encountered in previous perturbative 
calculations of l" 1 quarkonium decays are found to be explicitly cancelled, and a finite result for the 
decay width to order al is then obtained. 
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The study of heavy quarkonium plays an important role in the understanding of quantum chromody- 
namics (QCD). In recent years, since the E760 collaboration at Fermilab jl], |2| observed charmonium states 
via resonant pp-annihilation and measured their masses, widths, and branching fractions with unprecedented 
precision, problems about decay and production of heavy quarkonium have aroused much interest of people. 
In particular, with the observation of l" 1 charmonium state h c in experiment Jl], the study of its properties 
such as mass and decay width becomes very interesting both theoretically and experimentally. In the earlier 
treatment of heavy quarkonium annihilation, quarkonium is only taken as a color-singlet bound state of a 
heavy quark Q and its antiquark Q. Calculation of decay rate is based on the assumption that the annihila- 
tion of Q and Q is a short-distance process which, because of asymptotic freedom of QCD, can be computed 
in perturbation theory, and all nonperturbative effects could be factored into a constant: the wavefunction 
at the origin or its derivative at the origin. Using this factorization assumption to calculate the hadronic 
decay width of P-wave quarkonium 1^ state, the infrared divergence appears in the limit of small binding 
energy || and this divergence can not be factored into the derivative of the wavefunction. It is interpreted 
as a signal that the decay rate is sensitive to the nonperturbative effects which can not be contained in the 
wavefunction. This implies that the earlier factorization assumption fails because it is not rigorous within 
the framework of QCD. 

Recently, Bodwin, Braaten and Lepage || have developed a factorization formalism that allows 
systematic calculations of inclusive decay rates and production cross section of heavy quarkonium to any 
order in QCD coupling constant a s and to any order in v 2 , where v is the typical relative velocity of the heavy 
quark. The factorization formula is based on the use of the effective field theory NRQCD (Nonrelativistic 
QCD). Using this rigorous factorization formalism, the decay rate can be written as a sum of a set of long- 
distance nonperturbative matrix elements of which each is multiplied by a short-distance coefficient which 
can be calculated in perturbative QCD. This approach is successful in the study of many processes about 
decay and production of heavy quarkonium g . In this paper, we apply NRQCD to the hadronic decay of l" 1 
quarkonium and compute the QCD radiative corrections to its annihilation rate from both the color-singlet 
and color-octet QQ Fock states. We will show the explicit cancellation of the infrared divergence between 
the color-singlet and color-octet QQ components, and give a complete result at next to leading order in a s 
and leading order in v 2 . 

In NRQCD, the effect of annihilation can be taken into account by adding 4-fermion operators to 
NRQCD Lagrangian: 



where the sum is over all possible local 4-fermion operators O n that annihilate and create a QQ pair, and 




(1) 



n 
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d n is the scaling dimension of O n . The short distance coefficients f n (ot s ) can be computed by matching 
perturbative amplitudes for QQ scattering in NRQCD with the corresponding amplitudes in full QCD. The 
annihilation rate of a quarkonium state H to light hadrons (LH) can be written as 

T(H -» LH) = 21m < H\S£ 4 - fermion \H > (2) 

At a given order in v 2 , the number of matrix elements can be reduced to a finite number by using velocity 
scaling rules for the matrix elements Q ||. These scaling rules consist of that for the operators and for the 
probabilities of the Fock states that give the leading contributions to the matrix elements. 

At leading order in v 2 , the decay width of l" 1 quarkonium state can be written as 

r(l+- -> LH) = 2/m/i( 1 P 1 )iT 1 + 2Imf s CS )H s + 0{v 2 T), (3) 

where two nonperturbative parameters Hi and Hg can be defined rigorously in terms of matrix elements of 
a color-singlet and a color-octet 4-fermion operator in NRQCD 

< l+-|0 1 ( 1 P 1 )|l+-> 
Hl - m* ' 

< 1+-|Q 8 ( 1 5q)|1+- > 

#8 — 2 ' 

TTl 

where 

1 CP 1 ) = ^(- l -B)x-X + (-^m, 

O 8 ( 1 S )=iiJ + T a x -x + T a iP, 

where D is the space part of the covariant derivative D 11 and T a (a = 1, • • • , N 2 — l) is the SU (N c ) color matrix, 
and if) and x + are t ne fields with two components for quark Q and antiquark Q in NRQCD. Here including 
Hg is due to the fact that for decays of l" 1 heavy quarkonium, an S-wave color-octet QQ component in 
the wavefunction will contribute at the same order in v as the P-wave color-singlet QQ component, because 
the probability for annihilating an S-wave color-octet QQ state is proportional to v° while this state has a 
probability at order of v 2 to be transmitted into a P-wave color-singlet QQ state through the emission of 
soft gluon, whereas the dominate Fock state of l -1 quarkonium is \QQ > with QQ pair in a color-singlet 
1 P\ state, and the probability for annihilation through P-wave color singlet QQ is proportional to v 2 . 

In the following we first calculate the coefficients of non-perturbative matrix elements Hi and H^ 
to order by matching the imaginary part of perturbative scattering amplitude of QQ — > QQ in full QCD 
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with that in NRQCD, and then derive the formula of l -1 quarkonium decay width to an accuracy of next- 
to-leading order in a s . Using phenomenological parameters H± and Hg determined from other processes, we 
finally give an approximate numerical estimate of the decay width. 

We first calculate the imaginary part of QQ forward scattering amplitude ImM in full QCD. For 
convenience, we consider QQ scattering in the center of momentum frame with the the momenta of the heavy 
quarks and antiquarks small compared to the heavy quark mass. We take the incoming Q and Q to have 
momenta p and —p, while the outgoing Q and Q have momenta p 1 and —p 1 . By the conservation of energy, 
we have \p?\ = \p\ = p. In order to compare with the result in NRQCD, following [|], in the expression 
of ImM to be calculated in full perturbative QCD, we write the 4-component Dirac spinors in the Dirac 
representation in terms of 2-component Pauli spinors via the substitutions, 

where E = y/ m 2 + p 2 , £ and 77 are 2-component spinors with color indices suppressed. The Dirac spinors 
u{j?) and v(-p') have similar expressions in terms of Pauli spinors £' and rf . The spinors (Q) and (^) 
represent fermion states with standard nonrelativistic normalization. 



(a) 




(b) 



Fig.l Feynman diagrams contributing to C: 



■full QCD 

(T««r")(i®i) 



at order 



It is known from Q that to leading order in a s , only the coefficient /m/s( 1 5o) in (||) does not vanish 
and therefore only the color-octet matrix element Hg contributes to the decay width. In the S-wave case, 
we expand the annihilation amplitude ImM. in terms of velocities v = p/E and v' — p 1 jm only to leading 
order, and reduce ImM to four terms 



ImM 

1 «)T< 1 )((T i (»(T i 



(G) 



In order to determine (Imfsi 1 So))o, we only consider the coefficient C(r a ®T a )(i®i) °f ^ ne term £' + T a r]'r] + T a (;. 
At order a 2 s , only two diagrams shown in Fig. 1(a) and Fig. 1(b) contribute to this coefficient. After decom- 
posing the spinors and expanding them to leading order in v and v 1 , we obtain the coefficient of the term 
£+T a rfr)+T a £; in ImM 

= - 2)(d - 3)<&(2) ' (7) 

where the two massless particle phase space $(2) in d = 4 — 2e dimension is integrated to give 

*(2) = i(^r r(1 - e) 



87rMm 2y r(2-2e)' 
A simple calculation leads to the final expression 

r fullQCD _ tt(^c ~ 4)o£ ,47tm 4 r(l - e) 

^ ( r-®r.)(i 8 i) " ANcm 2 1 4^2- J r(2-2e) (1 ejU 2ej ' (8) 
where the dimensionless coupling constant is defined as 



While in NRQCD the QQ forward scattering amplitude can be reproduced by 4-fermion operators in 
the effective lagrangian, and the corresponding term £' + T a r]'r] + T a t; in ImM comes from operator 0%(}Sq) : 

of which the coefficient is IHlhL^sl, By comparing it with (||), we get 

(I-/s(%))o = ^ (^f^a - 0(1 - 2e). (9) 
We keep e in the above expression for convenience of later calculations. In the limit e — > 0, we get 

(Im/ 8 ( 1 5 ))o - 7r( ^~ 4) "«. ( 10 ) 

which has been given in 0. Here the subscript "0" in the coefficient means only the result at leading order 
in a s is taken, and the width can be written as 

r(l+" - LH) = %iN 2 * N ~ 4) <x 2 s H 8 + 0(a s T). (11) 

In order to obtain the next-to-leading order result, we must take account of the effects coming from 
both color-octet component and color-singlet component of the quarkonium. In the following, we calculate 
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the coefficients Imf\{ 1 P\) to leading order in a s and Imfsi 1 So) to next-to-leading order in a 8) and then 
give the complete formular for the hadronic decay width of l P\ to order of at leading order of v 2 . 
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(a) (b) (c) (d) 

Fig. 2 Feynman diagrams contributing to Cl^l^^^^ at leading order in a s 



Via the same procedure as above, we consider the imaginary part of QQ scattering amplitude, and 
calculate the coefficient Imf\Q-P\) by matching a perturbative calculation in full QCD with the correspond- 
ing perturbative calculation in NRQCD. The P-wave case requires an expansion of the annihilation amplitude 
ImM up to the first power of relative momenta p and p 1 . At this order ImM can be written as 

ImM 

= CV'-t7(i®i)(i®i)w" • v£,' + t]'ii + £ + CV'.^Kg,!)^.®^)^' • v£' + ar)' ■ 7/ + ct£ 

+ C v n v ^ im){a i 0aj) ^ + v' ■ ar/'r) + v-a£ + C t ,/*„«(i®i)( ff < (8<r 4)£' + tf ■ ar}'r] + nf ■ <r£ 

+ C 0l<T * 9T . )(191) if • $e + T a r]'r)+T a Z + C^ <TW){a i^ ■ vt'+aT a r,' ■ r,+aT a £ 

+ C v n vHTa ^ Ta)(M) C + v' ■ aT a V ' V +v-aT a ^ + C v , jvi(Ta0Ta)(ai0aj) e+v-cW a r]'r)+v' ■ aT a C (12) 

The determination of Jm/i( Pi) only requires calculating the coefficient C$i .fin^^n^u of the term v' ■ 
v£,' + r]'ri + £ t . In full QCD to leading order in a s , only the diagrams in Fig. 2 contribute to this coefficient. 
After making a nonrelativistic expansion for ImM to first order in v and v', we get 

r full QCD 

(N*-4)C F g e d-3 {[( _ 320 + 96rf) 1 +(768 _ 2 40d + 4d 2 )( 1 ■ 1 ^ 



8iV 2 m 4 48(d - 1) 
1 



• ryi3 -,2 ,-y>2 3 ' 

12 1 2 



(-512 + 176rf - 4^)(-3— + -^3) + (64 - 32d)(-^ + i) 



1 



1 



1 



(-1344 + 404d - 2d z )^-^ + (624 - 208rf + 2d z ){^— + ? ) 

- x^a^ X1X2 



l J '2 

(-64 + 32d)(i + + (-168 + 68d - 3d 2 )—!— ] 

X^ ^2 £1X2 

(two other permutations)} d$> (3). 



(13) 
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Here Xi — ki/m (i = 1,2,3) and ki denote the energies of the final-state gluons. The massless three-body 
phase space can be written as 

After performing the integration for two invariant X\ and X2, we obtain 

r fuii Q CD _ {N 2 -4)C F a 3 4tt 2 2 3 1 - 2e 1 7tt 2 - 94 

- 3iV 2 m 2 l 4m 2j ^ J T(2 - 2e) 1 2e /fl 48 ' 

{Imh( l S )) 4C F a S l I 1 4tt^ 7^-118 

where 7^ = 0.577 is the Euler constant. Comparing (|l4|) with the result obtained in ref. || which is regularized 
by the binding energy of QQ pair, we find that if making the substitution In^j — ► — , the two results 
have the same divergent terms, but their finite terms are different due to different regularization schemes. 
Here e is the binding energy of QQ pair, which is defined as 

e _ Am 2 - M 2 
m 4to 2 ' 

where M is the mass of QQ bound state. Here we control the infrared divergence using on-shell dimensional 
regularization, because off-shell binding energy regularization scheme will break manifest gauge invariance 
and conventional treatment of NRQCD is exact only for on-shell amplitudes. However, after taking account 
of the contribution from color-octet QQ component we will find that the coefficient /to/i( 1 Pi) is infrared 
finite and the final result is independent of infrared regularization scheme. 




Fig. 3 Feynman diagram contributing to Cv.^n^nnsii) through the operator Oi( 1 Pi) 



G 





Fig. 4 Feynman diagrams contributing to C^^n^^n®!) through the operator Og( 1 5o) 



In NRQCD, the QQ forward scattering amplitude can be reproduced by operators in bL^-f erm i on . 
When working at order q^, there are two 4-fermion operators which contribute to the coefficient C$i .^i^ijn&u 
of the term v' ■ v£' + r]'r] + !; in ImM, which are 

5£.4-fermion = ~ ^ Pi) + ^ - °^ Q$ ( 1 5'o) ■ 



The color singlet operator Oi( 1 Pi) contributes through the tree diagram in Fig. 3 which contains a 4-fermion 
vertex corresponding to Oi( 1 Pi) , and the result is 



ImM F iq.3 = 



-z V ■ Vt, T] T) f. 



(15) 



Since Imfg( So) is already known to be of order a 2 s , it is necessary to compute the contribution of the 
operator Og( 1 S'o) to an accuracy of a s . It is obvious that this contribution only comes from one-loop 
diagrams in Fig.4(a)-(d) which contain a 4-fermion vertex corresponding to Og( So), and these one-loop 
figures cause the transition from a color octet QQ into a color singlet QQ. The overall contribution of 
diagrams in Fig 4 is 



ImM.Fig.4 



Imf 8 ( 1 S )4:C F a s , 11 ^iAr 

2 77J 1E+ In f 

m z 3iV c 7r 2 e/_R 4m z 



2 ec/v 



3iV c 7r 1 2 y e IR 



■) 



(16) 



where is the IR (infrared) divergence and jim is the corresponding scale, while is the UV (ultraviolet) 

divergence and /ijjv is the corresponding scale. After the renormalization of operator Os^So) in the MS 
scheme the result is free from UV divergence, but the IR divergence still remains and it represents the 
nonperturbative nature of the annihilation amplitude. To order o? s , Imfs^So) on the right hand side of ( |l6| ) 
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must be taken as (Jm/s( 1 S'o))o, and then we obtain 

n NRQCD _ Imf l ( 1 P l ) (Ito/ 8 ( 1 5'o)) AC F a s r 11 ^Ar \ , , , 17 n 

From (0) and (0), we find that the coefficients of IR divergence are the same. It is clear that the IR 
divergence appearing in ([l4]) is proportional to the probability of transition between a color-singlet QQ pair 
and a color-octet QQ pair by the emission of soft gluon. This is the nonperturbative effect and must be 
factored into the long-distance matrix elements which have been defined explicitly in NRQCD. Comparing 
(|l4j) with ( |l7| ) and using (^), the finite coefficient Im/i( 1 Pi) is found to be 

Obviously the previously encountered IR divergence has been canceled and factored into the nonperturbative 
matrix element. The operator Os^Sq) satisfies the evolution equation 

M 5 = — 2°n Pi), (19) 

which has been derived in Q . We have neglected the subscript "£/ V" in /u and we will keep this notation in 
our work. 
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Fig. 5 Representive diagrams contributing to the first-order radiative correction to 

r ,full QCD 

U (T»ig)T°)(l(81) 



We have derived the coefficient Imfg^So) to leading order in a s . In order to get the result to 
next-to-leading order, we must consider the imaginary part of scattering amplitude of QQ pair to order in 
in full QCD. The diagrams which contribute to the coefficient of the term £' + T a 7/?/ + T a £ in ImM. to 
next-to-leading order in a s are shown in Fig. 5. We only give the representative diagrams and neglect the 
diagrams which give the same result as some of those in Fig. 5. The contribution from each diagram in terms 
of the unrenormalized coupling constant has in general the following form: 



(ImfsC S ))o a s 



f(e)A(" diagram") 



(20) 
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Table 1: C^aJ^^iv,!) from individual diagrams shown in Fig. 5 





Contribution from 


Contribution from 


Diagram 


two-particle cut 


three-particle cut 


(a) 


(C7 F -^)(fl + i i-+2/n2-2) 





(b) 


C F (-^-+3ln2-l) 





(c) 







(d) 







(e) 


C a[^7 - + + 2 - /n2 + fl] 

L ZtJ7V Z s t 6 ' _LZ J 


^A[h{-k + h)iR + h — if] 


(f) 






(g) 


"/ + 2 V ln m ' 

3errv r 3 lu 2m 


2 V In mi 8n/ 
3^! m 2m 9 


(h) 








(i) 








0) 


%[-(^ + 7)/«-2 + 2Zn2+^f] 


+ + 9 


(k) + (I) 





2 ^e/jf 3 ' 


(m) + (n) 









The values are normalized to {Imh £ 2 So))o /(e) ^. 



Here 6 = (4 - d)/2, /(e) = (^) e r(l + e), <7 F = C A = 7V C , 

n / stands for the number of light flavors 



with 

The imaginary part of these diagrams receives contributions from two-gluon cut, three-gluon cut and a 
"light" quark-antiquark pair plus one-gluon cut. The contribution of each individual diagram is calculated 
in Feynman gauge. Hence we have to add a ghost contribution both to two-gluon cut and to three-gluon cut 
in the diagram of Fig.5f. Our results for the contributions from individual diagrams are listed in table 1. 

Divergences show up in the intermediate steps of the calculation, the dimensional regularization 
procedure is used by going to d dimensions and introducing a scale /U through the standard replacement 

of the bare coupling constant g — > g^ d ~ 4 ^ 2 . Manifest gauge invariance and massless particle kinematics 
greatly simplify the calculations. The origin of the e = poles is specified in the table by the subscripts 
UV and IR. In the table we give the regularized und nrenormalized results for these diagrams, which show 
a l/(d — 4) divergence and a finite part. 

The overall result for the unrenormalizcd first-order radiative correction to the coefficient C(T"®T a )(i®i) 
in full QCD can be obtained by summing up all different individual contributions, and reads 

r full QCD 

U (T»ig>T»)(lig>l) 
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(Jm/8 ? 0)) ° {l + ^f{eW- + (C F - £±)£ + A]}, (21) 



where 



/!,i = ( ^ )£r(1 + £) ' 



60 = ^(11^-271/), 

A = C F {— 5) + ^(—0 - _) - ~n f , 

with Cf = -%y — > Ca — N c . The term with - indicates the Coulomb singularity which arises from the 
Coulomb exchange of the gluon between quark and antiquark in Fig. 5(a). It is sensitive to the long-distance 
nonperturbative effect. 

We find that the cancellation of the infrared divergences occurs in the overall result in spite of 
the fact that the individual cuts do not. This is the same as the corresponding color singlet coefficient 
^(i®iKi®i)' ^ s a ma tt er of fact it may be interesting to know that the infrared divergences of Fig.5a,c,k,l 
cancel each other, which is different from the color singlet coefficient, where the cancellation occurs between 
Fig. 5a and Fig. 5c, as well as between Fig. 5k, Fig. 51, Fig. 5m and Fig.5n respectively J|, ||. 

Now we renormalize the coupling constant in the MS scheme with 

" "- lfS (i_<!6 (i + z n4 7r-7 E )), 



and find 



7T 7T 7T 6 



r full QCD 



^i) oHl + ^ ((Cr _^ + 4W „|_ + , 1} , <*> 



where we have suppressed the superscript MS in a s 



to thank Dr. Petrelli et al. to point out the problem of this term for our old version in their recent paper 
ticp-ph/9707223 Now in this new version it is important to note that the two independent calculations are found to give a 



identical result. 
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Fig. 6 Feynman diagrams contributing to C(T<*®T a )(i®i) to next-to-leading order of a s in 

NRQCD 



In order to determine Imf^ 1 So), we must calculate the corresponding contribution of 5C^-f e 



to C, 



full QCD 
(T"®T")(1®1) 



jn in NRQCD to next-to-leading order in a s . The relevant Feynman diagrams are shown 



in Fig. 6. They contain a four-fermion vertex that corresponds to the term ip + T a xx + T a ip in the effective 
Lagrangian. In the limit v — > 0, only Fig. 6(b) and Fig. 6(c), which include Coulomb exchange of the gluon, 
contribute at next-to-leading order. Fig. 6(a) gives leading order result 



ImM 



6(a) 



(23) 



The contribution from Fig. 6(b) is 



ImM 



Imfs^Sp) C A ira s i 1 /+ , + 



where the imaginary part arises because the incoming quark and antiquark can scatter on shell before being 
annihilated by the 4-fermion operator. The contribution from Fig. 6(c) is 

ImM 6{b) = !^A {Cf - + £( J_ - 1E + lnbr - 2ln™rW{. (25) 

m z 2 4u 7r e/fl, 



Add the contributions from Fig. 6(a), (b),(c) together, we obtain the complete result for C(ti®t°-)(i®i) to 
order of a 3 s in NRQCD 



r NRQCD 



Jm/ 8 ( 1 S'o) r a s c a,k 2 

2 1 H \ C F TTITT ■ 



(26) 



Comparing ( |22j ) and (Eq), we can read off the imaginary part of /s( 1 <S'o) to next-to-leading order in 



7m/ 8 ( 1 5 ) - 



(27) 
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Note that the factorization approach reproduces the standard prescription of simply dropping the l/v terms 
in the perturbatively calculated annihilation rate. It is clear that the Coulomb singularity can be factored 
into the nonperturbative part trivially in this factorization formula. 

Having derived the coefficients Imfi^Pi) and /m/g( 1 So), we finally come to the overall result for 
hadronic decay of l" 1 quarkonium state to next-to-leading order in a s at leading order of v 2 . Substituting 
@ and (§7|) into (§), we get 



TV1+- rm 2(iV c 2 -4)C F a^ 7^-118 „ 



+ ^^(l + ^Wn^i))%). (28) 
2l\ c tt 2m 

Working to all orders in a s (/i), the final result is independent of /i, since the coefficients depend on fj, in such 
a way that they will cancel the [i dependence of matrix elements. 

Now we apply our above result to the charmonium system to study the decay width of h c . In ( p8| ) 
making a choice of [i = m c and taking N c = 3, rif = 3 we obtain 

Cy ( 777 I 

T(h c -> LH) = -0.16al(m c )Hi + 2.62a 2 Jm c )(l + 7.10 " K cJ )H 8 (29) 

7T 

It is interesting to note that the contribution of color singlet component is negative and the QCD radiative 
correction from color octet component is very large. The matrix elements Hi and Hg have been defined 
explicitly in NRQCD and they are difficult to derive from first principles of QCD. People have tried to 



compute them using lattice simulations 10 . In practice they can be determined phenomenologically. Heavy 
quark spin symmetry provides approximate relations between them and the corresponding two parameters 
in the expressions of decay widths of P-wave triplet Xcj(J = 0, 1, 2) states. At leading order of v 2 , they are 
equal respectively. A rough estimate of Hi and H% have been given in ]l2| by comparing the theoretical 
result of Xcj decay to order with experimental data. There they don't give the coefficients of H$ at order 
a^, because in Xcj decays, the contributions of color-octet component are the same for J=0,l,2, and can be 
treated as just one parameter. However it is not Hg, mentioned above. Here, as an approximation, using the 
estimated value for Hi and H$ to the order of a 2 in 0, 

Hi = 15.3 ± 3.7Mev, H s = 3.26 ± 0.73Meu, 

a s {m c ) =0.25 ±0.02, 

we roughly get T = 0.80 ± 0.20Mev. A more reliable estimate will be obtained with a complete theoretical 



result for the \cj decay width to order a 
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In this work we use a general factorization formula which is based on NRQCD to calculate the 
annihilation rate of l" 1 quarkonium, we sec that the infrared divergence appearing in previous calculations 
can be factored into the long-distance perturbative matrix element rigorously. Our result is also free from 
the Coulomb singularity. The corresponding case of XcJ production through gluon fragmentation has been 
studied in JIl| . It is clear from our calculation that the failure of previous factorization assumption is due 
to the fact that only color singlet component was considered and all contributions from color octet were 
neglected. In that sense the previous result is incomplete and therefore the infrared divergence may appear 
in some cases such as the annihilation and production of P-wave quarkonium even at leading order in v 2 . 
Our calculation shows that the rigorous factorization formula can separate short-distance perturbative effects 
from long-distance nonperturbative effects correctly and can therefore provide a systematical calculation for 
quarkonium decay and production to any order in a s and in v 2 , because it is based on a solid theoretical 
foundation. 
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